home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / gsl-randist.c < prev    next >
Encoding:
C/C++ Source or Header  |  2001-08-18  |  10.0 KB  |  394 lines
  1. /* randist/gsl-randist.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdio.h>
  22. #include <stdlib.h>
  23. #include <math.h>
  24. #include <string.h>
  25.  
  26. #include <gsl/gsl_randist.h>
  27. #include <gsl/gsl_rng.h>
  28. #include <gsl/gsl_test.h>
  29.  
  30. void error (const char * s);
  31.  
  32.  
  33. int
  34. main (int argc, char *argv[])
  35. {
  36.   size_t i,j;
  37.   size_t n = 0;
  38.   double mu = 0, nu = 0, nu1 = 0, nu2 = 0, sigma = 0, a = 0, b = 0, c = 0;
  39.   double zeta = 0, sigmax = 0, sigmay = 0, rho = 0;
  40.   double p = 0;
  41.   double x = 0, y =0, z=0  ;
  42.   unsigned int N = 0, t = 0, n1 = 0, n2 = 0 ;
  43.   unsigned long int seed = 0 ;
  44.   const char * name ;
  45.   gsl_rng * r ;
  46.  
  47.   if (argc < 4) 
  48.     {
  49.       printf (
  50. "Usage: gsl-randist seed n DIST param1 param2 ...\n"
  51. "Generates n samples from the distribution DIST with parameters param1,\n"
  52. "param2, etc. Valid distributions are,\n"
  53. "\n"
  54. "  beta\n"
  55. "  binomial\n"
  56. "  bivariate-gaussian\n"
  57. "  cauchy\n"
  58. "  chisq\n"
  59. "  dir-2d\n"
  60. "  dir-3d\n"
  61. "  dir-nd\n"
  62. "  erlang\n"
  63. "  exponential\n"
  64. "  exppow\n"
  65. "  fdist\n"
  66. "  flat\n"
  67. "  gamma\n"
  68. "  gaussian-tail\n"
  69. "  gaussian\n"
  70. "  geometric\n"
  71. "  gumbel1\n"
  72. "  gumbel2\n"
  73. "  hypergeometric\n"
  74. "  laplace\n"
  75. "  landau\n"
  76. "  levy\n"
  77. "  levy-skew\n"
  78. "  logarithmic\n"
  79. "  logistic\n"
  80. "  lognormal\n"
  81. "  negative-binomial\n"
  82. "  pareto\n"
  83. "  pascal\n"
  84. "  poisson\n"
  85. "  rayleigh-tail\n"
  86. "  rayleigh\n"
  87. "  tdist\n"
  88. "  ugaussian-tail\n"
  89. "  ugaussian\n"
  90. "  weibull\n") ;
  91.       exit (0);
  92.     }
  93.  
  94.   argv++ ; seed = atol (argv[0]); argc-- ;
  95.   argv++ ; n = atol (argv[0]); argc-- ;
  96.   argv++ ; name = argv[0] ; argc-- ; argc-- ;
  97.  
  98.   gsl_rng_env_setup() ;
  99.  
  100.   if (gsl_rng_default_seed != 0) {
  101.     fprintf(stderr, 
  102.         "overriding GSL_RNG_SEED with command line value, seed = %ld\n", 
  103.         seed) ;
  104.   }
  105.   
  106.   gsl_rng_default_seed = seed ;
  107.  
  108.   r = gsl_rng_alloc(gsl_rng_default) ;
  109.  
  110.  
  111. #define NAME(x) !strcmp(name,(x))
  112. #define OUTPUT(x) for (i = 0; i < n; i++) { printf("%g\n", (x)) ; }
  113. #define OUTPUT1(a,x) for(i = 0; i < n; i++) { a ; printf("%g\n", x) ; }
  114. #define OUTPUT2(a,x,y) for(i = 0; i < n; i++) { a ; printf("%g %g\n", x, y) ; }
  115. #define OUTPUT3(a,x,y,z) for(i = 0; i < n; i++) { a ; printf("%g %g %g\n", x, y, z) ; }
  116. #define INT_OUTPUT(x) for (i = 0; i < n; i++) { printf("%d\n", (x)) ; }
  117. #define ARGS(x,y) if (argc != x) error(y) ;
  118. #define DBL_ARG(x) if (argc) { x=atof((++argv)[0]);argc--;} else {error( #x);};
  119. #define INT_ARG(x) if (argc) { x=atoi((++argv)[0]);argc--;} else {error( #x);};
  120.  
  121.   if (NAME("bernoulli"))
  122.     {
  123.       ARGS(1, "p = probability of success");
  124.       DBL_ARG(p)
  125.       INT_OUTPUT(gsl_ran_bernoulli (r, p));
  126.     }
  127.   else if (NAME("beta"))
  128.     {
  129.       ARGS(2, "a,b = shape parameters");
  130.       DBL_ARG(a)
  131.       DBL_ARG(b)
  132.       OUTPUT(gsl_ran_beta (r, a, b));
  133.     }
  134.   else if (NAME("binomial"))
  135.     {
  136.       ARGS(2, "p = probability, N = number of trials");
  137.       DBL_ARG(p)
  138.       INT_ARG(N)
  139.       INT_OUTPUT(gsl_ran_binomial (r, p, N));
  140.     }
  141.   else if (NAME("cauchy"))
  142.     {
  143.       ARGS(1, "a = scale parameter");
  144.       DBL_ARG(a)
  145.       OUTPUT(gsl_ran_cauchy (r, a));
  146.     }
  147.   else if (NAME("chisq"))
  148.     {
  149.       ARGS(1, "nu = degrees of freedom");
  150.       DBL_ARG(nu)
  151.       OUTPUT(gsl_ran_chisq (r, nu));
  152.     }
  153.   else if (NAME("erlang"))
  154.     {
  155.       ARGS(2, "a = scale parameter, b = order");
  156.       DBL_ARG(a)
  157.       DBL_ARG(b)
  158.       OUTPUT(gsl_ran_erlang (r, a, b));
  159.     }
  160.   else if (NAME("exponential"))
  161.     {
  162.       ARGS(1, "mu = mean value");
  163.       DBL_ARG(mu) ;
  164.       OUTPUT(gsl_ran_exponential (r, mu));
  165.     }
  166.   else if (NAME("exppow"))
  167.     {
  168.       ARGS(2, "a = scale parameter, b = power (1=exponential, 2=gaussian)");
  169.       DBL_ARG(a) ;
  170.       DBL_ARG(b) ;
  171.       OUTPUT(gsl_ran_exppow (r, a, b));
  172.     }
  173.   else if (NAME("fdist"))
  174.     {
  175.       ARGS(2, "nu1, nu2 = degrees of freedom parameters");
  176.       DBL_ARG(nu1) ;
  177.       DBL_ARG(nu2) ;
  178.       OUTPUT(gsl_ran_fdist (r, nu1, nu2));
  179.     }
  180.   else if (NAME("flat"))
  181.     {
  182.       ARGS(2, "a = lower limit, b = upper limit");
  183.       DBL_ARG(a) ;
  184.       DBL_ARG(b) ;
  185.       OUTPUT(gsl_ran_flat (r, a, b));
  186.     }
  187.   else if (NAME("gamma"))
  188.     {
  189.       ARGS(2, "a = order, b = scale");
  190.       DBL_ARG(a) ;
  191.       DBL_ARG(b) ;
  192.       OUTPUT(gsl_ran_gamma (r, a, b));
  193.     }
  194.   else if (NAME("gaussian"))
  195.     {
  196.       ARGS(1, "sigma = standard deviation");
  197.       DBL_ARG(sigma) ;
  198.       OUTPUT(gsl_ran_gaussian (r, sigma));
  199.     }
  200.   else if (NAME("gaussian-tail"))
  201.     {
  202.       ARGS(2, "a = lower limit, sigma = standard deviation");
  203.       DBL_ARG(a) ;
  204.       DBL_ARG(sigma) ;
  205.       OUTPUT(gsl_ran_gaussian_tail (r, a, sigma));
  206.     }
  207.   else if (NAME("ugaussian"))
  208.     {
  209.       ARGS(0, "unit gaussian, no parameters required");
  210.       OUTPUT(gsl_ran_ugaussian (r));
  211.     }
  212.   else if (NAME("ugaussian-tail"))
  213.     {
  214.       ARGS(1, "a = lower limit");
  215.       DBL_ARG(a) ;
  216.       OUTPUT(gsl_ran_ugaussian_tail (r, a));
  217.     }
  218.   else if (NAME("bivariate-gaussian"))
  219.     {
  220.       ARGS(3, "sigmax = x std.dev., sigmay = y std.dev., rho = correlation");
  221.       DBL_ARG(sigmax) ;
  222.       DBL_ARG(sigmay) ;
  223.       DBL_ARG(rho) ;
  224.       OUTPUT2(gsl_ran_bivariate_gaussian (r, sigmax, sigmay, rho, &x, &y), 
  225.           x, y);
  226.     }
  227.   else if (NAME("dir-2d"))
  228.     {
  229.       OUTPUT2(gsl_ran_dir_2d (r, &x, &y), x, y);
  230.     }
  231.   else if (NAME("dir-3d"))
  232.     {
  233.       OUTPUT3(gsl_ran_dir_3d (r, &x, &y, &z), x, y, z);
  234.     }
  235.   else if (NAME("dir-nd"))
  236.     {
  237.       double *xarr;  
  238.       ARGS(1, "n1 = number of dimensions of hypersphere"); 
  239.       INT_ARG(n1) ;
  240.       xarr = (double *)malloc(n1*sizeof(double));
  241.  
  242.       for(i = 0; i < n; i++) { 
  243.         gsl_ran_dir_nd (r, n1, xarr) ; 
  244.         for (j = 0; j < n1; j++) { 
  245.           if (j) putchar(' '); 
  246.           printf("%g", xarr[j]) ; 
  247.         } 
  248.         putchar('\n'); 
  249.       } ;
  250.  
  251.       free(xarr);
  252.     }  
  253.   else if (NAME("geometric"))
  254.     {
  255.       ARGS(1, "p = bernoulli trial probability of success");
  256.       DBL_ARG(p) ;
  257.       INT_OUTPUT(gsl_ran_geometric (r, p));
  258.     }
  259.   else if (NAME("gumbel1"))
  260.     {
  261.       ARGS(2, "a = order, b = scale parameter");
  262.       DBL_ARG(a) ;
  263.       DBL_ARG(b) ;
  264.       OUTPUT(gsl_ran_gumbel1 (r, a, b));
  265.     }
  266.   else if (NAME("gumbel2"))
  267.     {
  268.       ARGS(2, "a = order, b = scale parameter");
  269.       DBL_ARG(a) ;
  270.       DBL_ARG(b) ;
  271.       OUTPUT(gsl_ran_gumbel2 (r, a, b));
  272.     }
  273.   else if (NAME("hypergeometric"))
  274.     {
  275.       ARGS(3, "n1 = tagged population, n2 = untagged population, t = number of trials");
  276.       INT_ARG(n1) ;
  277.       INT_ARG(n2) ;
  278.       INT_ARG(t) ;
  279.       INT_OUTPUT(gsl_ran_hypergeometric (r, n1, n2, t));
  280.     }
  281.   else if (NAME("laplace"))
  282.     {
  283.       ARGS(1, "a = scale parameter");
  284.       DBL_ARG(a) ;
  285.       OUTPUT(gsl_ran_laplace (r, a));
  286.     }
  287.   else if (NAME("landau"))
  288.     {
  289.       ARGS(0, "no arguments required");
  290.       OUTPUT(gsl_ran_landau (r));
  291.     }
  292.   else if (NAME("levy"))
  293.     {
  294.       ARGS(2, "c = scale, a = power (1=cauchy, 2=gaussian)");
  295.       DBL_ARG(c) ;
  296.       DBL_ARG(a) ;
  297.       OUTPUT(gsl_ran_levy (r, c, a));
  298.     }
  299.   else if (NAME("levy-skew"))
  300.     {
  301.       ARGS(3, "c = scale, a = power (1=cauchy, 2=gaussian), b = skew");
  302.       DBL_ARG(c) ;
  303.       DBL_ARG(a) ;
  304.       DBL_ARG(b) ;
  305.       OUTPUT(gsl_ran_levy_skew (r, c, a, b));
  306.     }
  307.   else if (NAME("logarithmic"))
  308.     {
  309.       ARGS(1, "p = probability");
  310.       DBL_ARG(p) ;
  311.       INT_OUTPUT(gsl_ran_logarithmic (r, p));
  312.     }
  313.   else if (NAME("logistic"))
  314.     {
  315.       ARGS(1, "a = scale parameter");
  316.       DBL_ARG(a) ;
  317.       OUTPUT(gsl_ran_logistic (r, a));
  318.     }
  319.   else if (NAME("lognormal"))
  320.     {
  321.       ARGS(2, "zeta = location parameter, sigma = scale parameter");
  322.       DBL_ARG(zeta) ;
  323.       DBL_ARG(sigma) ;
  324.       OUTPUT(gsl_ran_lognormal (r, zeta, sigma));
  325.     }
  326.   else if (NAME("negative-binomial"))
  327.     {
  328.       ARGS(2, "p = probability, a = order");
  329.       DBL_ARG(p) ;
  330.       DBL_ARG(a) ;
  331.       INT_OUTPUT(gsl_ran_negative_binomial (r, p, a));
  332.     }
  333.   else if (NAME("pareto"))
  334.     {
  335.       ARGS(2, "a = power, b = scale parameter");
  336.       DBL_ARG(a) ;
  337.       DBL_ARG(b) ;
  338.       OUTPUT(gsl_ran_pareto (r, a, b));
  339.     }
  340.   else if (NAME("pascal"))
  341.     {
  342.       ARGS(2, "p = probability, n = order (integer)");
  343.       DBL_ARG(p) ;
  344.       INT_ARG(N) ;
  345.       INT_OUTPUT(gsl_ran_pascal (r, p, N));
  346.     }
  347.   else if (NAME("poisson"))
  348.     {
  349.       ARGS(1, "mu = scale parameter");
  350.       DBL_ARG(mu) ;
  351.       INT_OUTPUT(gsl_ran_poisson (r, mu));
  352.     }
  353.   else if (NAME("rayleigh"))
  354.     {
  355.       ARGS(1, "sigma = scale parameter");
  356.       DBL_ARG(sigma) ;
  357.       OUTPUT(gsl_ran_rayleigh (r, sigma));
  358.     }
  359.   else if (NAME("rayleigh-tail"))
  360.     {
  361.       ARGS(2, "a = lower limit, sigma = scale parameter");
  362.       DBL_ARG(a) ;
  363.       DBL_ARG(sigma) ;
  364.       OUTPUT(gsl_ran_rayleigh_tail (r, a, sigma));
  365.     }
  366.   else if (NAME("tdist"))
  367.     {
  368.       ARGS(1, "nu = degrees of freedom");
  369.       DBL_ARG(nu) ;
  370.       OUTPUT(gsl_ran_tdist (r, nu));
  371.     }
  372.   else if (NAME("weibull"))
  373.     {
  374.       ARGS(2, "a = scale parameter, b = exponent");
  375.       DBL_ARG(a) ;
  376.       DBL_ARG(b) ;
  377.       OUTPUT(gsl_ran_weibull (r, a, b));
  378.     }
  379.   else
  380.     {
  381.       fprintf(stderr,"Error: unrecognized distribution: %s\n", name) ;
  382.     }
  383.  
  384.   return 0 ;
  385. }
  386.  
  387.  
  388. void
  389. error (const char * s)
  390. {
  391.   fprintf(stderr, "Error: arguments should be %s\n",s) ;
  392.   exit (EXIT_FAILURE) ;
  393. }
  394.